Stokes I Simultaneous Image and Instrument Modeling
In this tutorial, we will create a preliminary reconstruction of the 2017 M87 data on April 6 by simultaneously creating an image and model for the instrument. By instrument model, we mean something akin to self-calibration in traditional VLBI imaging terminology. However, unlike traditional self-cal, we will at each point in our parameter space effectively explore the possible self-cal solutions. This will allow us to constrain and marginalize over the instrument effects, such as time variable gains.
To get started we load Comrade.
using Comrade Activating project at `~/work/Comrade.jl/Comrade.jl/examples`using Pyehtim
using LinearAlgebraFor reproducibility we use a stable random number genreator
using StableRNGs
rng = StableRNG(42)StableRNGs.LehmerRNG(state=0x00000000000000000000000000000055)Load the Data
To download the data visit https://doi.org/10.25739/g85n-f134 First we will load our data:
obs = ehtim.obsdata.load_uvfits(joinpath(dirname(pathof(Comrade)), "..", "examples", "SR1_M87_2017_096_hi_hops_netcal_StokesI.uvfits"))Python: <ehtim.obsdata.Obsdata object at 0x7fc68c5e10f0>Now we do some minor preprocessing:
- Scan average the data since the data have been preprocessed so that the gain phases coherent.
- Add 1% systematic noise to deal with calibration issues that cause 1% non-closing errors.
obs = scan_average(obs.add_fractional_noise(0.01))Python: <ehtim.obsdata.Obsdata object at 0x7fc6d4f2ce80>Now we extract our complex visibilities.
dvis = extract_table(obs, ComplexVisibilities())EHTObservation{Float64,Comrade.EHTVisibilityDatum{Float64}, ...}
source: M87
mjd: 57849
frequency: 2.29070703125e11
bandwidth: 1.856e9
stations: [:AA, :AP, :AZ, :JC, :LM, :PV, :SM]
nsamples: 284
##Building the Model/Posterior
Now, we must build our intensity/visibility model. That is, the model that takes in a named tuple of parameters and perhaps some metadata required to construct the model. For our model, we will use a raster or ContinuousImage for our image model. The model is given below:
function sky(θ, metadata)
(;fg, c, σimg) = θ
(;ftot, K, meanpr, grid, cache) = metadata
# Transform to the log-ratio pixel fluxes
cp = meanpr .+ σimg.*c.params
# Transform to image space
rast = (ftot*(1-fg))*K(to_simplex(CenteredLR(), cp))
img = IntensityMap(rast, grid)
m = ContinuousImage(img, cache)
# Add a large-scale gaussian to deal with the over-resolved mas flux
g = modify(Gaussian(), Stretch(μas2rad(250.0), μas2rad(250.0)), Renormalize(ftot*fg))
return m + g
endsky (generic function with 1 method)Unlike other imaging examples (e.g., Imaging a Black Hole using only Closure Quantities) we also need to include a model for the instrument, i.e., gains as well. The gains will be broken into two components
- Gain amplitudes which are typically known to 10-20%, except for LMT, which has amplitudes closer to 50-100%.
- Gain phases which are more difficult to constrain and can shift rapidly.
function instrument(θ, metadata)
(; lgamp, gphase) = θ
(; gcache, gcachep) = metadata
# Now form our instrument model
gvis = exp.(lgamp)
gphase = exp.(1im.*gphase)
jgamp = jonesStokes(gvis, gcache)
jgphase = jonesStokes(gphase, gcachep)
return JonesModel(jgamp*jgphase)
endinstrument (generic function with 1 method)The model construction is very similar to Imaging a Black Hole using only Closure Quantities, except we include a large scale gaussian since we want to model the zero baselines. For more information about the image model please read the closure-only example. Let's discuss the instrument model Comrade.JonesModel. Thanks to the EHT pre-calibration, the gains are stable over scans. Therefore, we can model the gains on a scan-by-scan basis. To form the instrument model, we need our
- Our (log) gain amplitudes and phases are given below by
lgampandgphase - Our function or cache that maps the gains from a list to the stations they impact
gcache. - The set of
Comrade.JonesPairsproduced byjonesStokes
These three ingredients then specify our instrument model. The instrument model can then be combined with our image model cimg to form the total JonesModel.
Now, let's set up our image model. The EHT's nominal resolution is 20-25 μas. Additionally, the EHT is not very sensitive to a larger field of view. Typically 60-80 μas is enough to describe the compact flux of M87. Given this, we only need to use a small number of pixels to describe our image.
npix = 32
fovx = μas2rad(150.0)
fovy = μas2rad(150.0)7.27220521664304e-10Now let's form our cache's. First, we have our usual image cache which is needed to numerically compute the visibilities.
grid = imagepixels(fovx, fovy, npix, npix)
buffer = IntensityMap(zeros(npix, npix), grid)
cache = create_cache(NFFTAlg(dvis), buffer, BSplinePulse{3}())VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [1.2065151329523807e9 819926.8691406249 … 4.472903950222221e9 -15688.352240668402; -3.8185738483809524e9 -1.7123166484375e6 … -6.2063795222222224e7 119336.55902777778]), NFFTPlan with 284 sampling points for an input array of size(32, 32) and an output array of size(284,) with dims 1:2, ComplexF64[0.9304270810105981 - 0.1755525366839422im, 0.9999999857224657 - 6.37119066797087e-5im, 0.6968566352375541 - 0.5221791050707502im, 0.9304860627318433 + 0.17550202040945986im, 0.8028821861479044 + 0.3950138433443315im, 0.696909133455242 + 0.522149977275377im, 0.9999999859666728 - 6.454986502401626e-5im, 0.690100908057055 - 0.5259755617707915im, 0.9317353097306306 - 0.1646343527093735im, 0.9317937328057632 + 0.1645830409592267im … 0.8346090489457091 + 0.1795789988845693im, 0.834610800886739 + 0.17958571977595708im, 0.9181402138641586 + 0.09206131734480663im, 0.9120864418553256 + 0.2942607208902332im, 0.9120838428437319 + 0.29426742718619836im, 0.9559925804493471 + 0.09847981391058878im, 0.8881131676530709 + 0.289296040490046im, 0.9924358388500316 + 0.0028007013218525043im, 0.8881114735852544 + 0.28930284794286826im, 0.9999999999233905 + 7.3999334332420275e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])Second, we now construct our instrument model cache. This tells us how to map from the gains to the model visibilities. However, to construct this map, we also need to specify the observation segmentation over which we expect the gains to change. This is specified in the second argument to jonescache, and currently, there are two options
FixedSeg(val): Fixes the corruption to the valuevalfor all time. This is usefule for reference stationsScanSeg(): which forces the corruptions to only change from scan-to-scanTrackSeg(): which forces the corruptions to be constant over a night's observation
For this work, we use the scan segmentation for the gain amplitudes since that is roughly the timescale we expect them to vary. For the phases we use a station specific scheme where we set AA to be fixed to unit gain because it will function as a reference station.
gcache = jonescache(dvis, ScanSeg())
gcachep = jonescache(dvis, ScanSeg(); autoref=SEFDReference((complex(1.0))))
using VLBIImagePriorsNow we need to specify our image prior. For this work we will use a Gaussian Markov Random field prior Since we are using a Gaussian Markov random field prior we need to first specify our mean image. This behaves somewhat similary to a entropy regularizer in that it will start with an initial guess for the image structure. For this tutorial we will use a a symmetric Gaussian with a FWHM of 60 μas
fwhmfac = 2*sqrt(2*log(2))
mpr = modify(Gaussian(), Stretch(μas2rad(50.0)./fwhmfac))
imgpr = intensitymap(mpr, grid)2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓ X ∈ 32-element LinRange{Float64,...}
→ Y ∈ 32-element LinRange{Float64,...}
And data, 32×32 NamedDimsArray(::Matrix{Float64}, (:X, :Y)):
(-3.52247e-10) … (3.29522e-10) (3.52247e-10)
(-3.52247e-10) 6.37537e-8 1.32434e-7 6.37537e-8
(-3.29522e-10) 1.32434e-7 2.751e-7 1.32434e-7
(-3.06796e-10) 2.62014e-7 5.44273e-7 2.62014e-7
(-2.84071e-10) 4.93724e-7 1.0256e-6 4.93724e-7
(-2.61345e-10) 8.86092e-7 … 1.84065e-6 8.86092e-7
(-2.38619e-10) 1.51463e-6 3.14629e-6 1.51463e-6
⋮ ⋱ ⋮
(2.15894e-10) 2.46586e-6 5.12225e-6 2.46586e-6
(2.38619e-10) 1.51463e-6 3.14629e-6 1.51463e-6
(2.61345e-10) 8.86092e-7 … 1.84065e-6 8.86092e-7
(2.84071e-10) 4.93724e-7 1.0256e-6 4.93724e-7
(3.06796e-10) 2.62014e-7 5.44273e-7 2.62014e-7
(3.29522e-10) 1.32434e-7 2.751e-7 1.32434e-7
(3.52247e-10) 6.37537e-8 1.32434e-7 6.37537e-8Now since we are actually modeling our image on the simplex we need to ensure that our mean image has unit flux
imgpr ./= flux(imgpr)2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓ X ∈ 32-element LinRange{Float64,...}
→ Y ∈ 32-element LinRange{Float64,...}
And data, 32×32 Matrix{Float64}:
(-3.52247e-10) … (3.29522e-10) (3.52247e-10)
(-3.52247e-10) 6.38049e-8 1.3254e-7 6.38049e-8
(-3.29522e-10) 1.3254e-7 2.75321e-7 1.3254e-7
(-3.06796e-10) 2.62224e-7 5.4471e-7 2.62224e-7
(-2.84071e-10) 4.94121e-7 1.02642e-6 4.94121e-7
(-2.61345e-10) 8.86803e-7 … 1.84213e-6 8.86803e-7
(-2.38619e-10) 1.51585e-6 3.14882e-6 1.51585e-6
⋮ ⋱ ⋮
(2.15894e-10) 2.46784e-6 5.12636e-6 2.46784e-6
(2.38619e-10) 1.51585e-6 3.14882e-6 1.51585e-6
(2.61345e-10) 8.86803e-7 … 1.84213e-6 8.86803e-7
(2.84071e-10) 4.94121e-7 1.02642e-6 4.94121e-7
(3.06796e-10) 2.62224e-7 5.4471e-7 2.62224e-7
(3.29522e-10) 1.3254e-7 2.75321e-7 1.3254e-7
(3.52247e-10) 6.38049e-8 1.3254e-7 6.38049e-8and since our prior is not on the simplex we need to convert it to unconstrained or real space.
meanpr = to_real(CenteredLR(), Comrade.baseimage(imgpr))32×32 Matrix{Float64}:
-7.55422 -6.82317 -6.14085 -5.50727 … -6.14085 -6.82317 -7.55422
-6.82317 -6.09211 -5.4098 -4.77622 -5.4098 -6.09211 -6.82317
-6.14085 -5.4098 -4.72748 -4.0939 -4.72748 -5.4098 -6.14085
-5.50727 -4.77622 -4.0939 -3.46032 -4.0939 -4.77622 -5.50727
-4.92243 -4.19137 -3.50906 -2.87548 -3.50906 -4.19137 -4.92243
-4.38632 -3.65527 -2.97295 -2.33937 … -2.97295 -3.65527 -4.38632
-3.89895 -3.1679 -2.48558 -1.852 -2.48558 -3.1679 -3.89895
-3.46032 -2.72927 -2.04695 -1.41337 -2.04695 -2.72927 -3.46032
-3.07043 -2.33937 -1.65705 -1.02348 -1.65705 -2.33937 -3.07043
-2.72927 -1.99821 -1.3159 -0.682317 -1.3159 -1.99821 -2.72927
⋮ ⋱ ⋮
-3.07043 -2.33937 -1.65705 -1.02348 -1.65705 -2.33937 -3.07043
-3.46032 -2.72927 -2.04695 -1.41337 -2.04695 -2.72927 -3.46032
-3.89895 -3.1679 -2.48558 -1.852 … -2.48558 -3.1679 -3.89895
-4.38632 -3.65527 -2.97295 -2.33937 -2.97295 -3.65527 -4.38632
-4.92243 -4.19137 -3.50906 -2.87548 -3.50906 -4.19137 -4.92243
-5.50727 -4.77622 -4.0939 -3.46032 -4.0939 -4.77622 -5.50727
-6.14085 -5.4098 -4.72748 -4.0939 -4.72748 -5.4098 -6.14085
-6.82317 -6.09211 -5.4098 -4.77622 … -5.4098 -6.09211 -6.82317
-7.55422 -6.82317 -6.14085 -5.50727 -6.14085 -6.82317 -7.55422Now we can form our metadata we need to fully define our model.
metadata = (;ftot=1.1, K=CenterImage(imgpr), meanpr, grid, cache, gcache, gcachep)(ftot = 1.1, K = VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}([0.9944957386363662 -0.005326704545454548 … 0.005326704545454572 0.005504261363636381; -0.005326704545454548 0.9948393969941349 … 0.005160603005865089 0.005326704545454565; … ; 0.005326704545454572 0.005160603005865089 … 0.9948393969941348 -0.005326704545454546; 0.005504261363636381 0.005326704545454565 … -0.005326704545454546 0.9944957386363638], (32, 32)), meanpr = [-7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; … ; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; -7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378], grid = GriddedKeys{(:X, :Y)}
X: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
Y: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
, cache = VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [1.2065151329523807e9 819926.8691406249 … 4.472903950222221e9 -15688.352240668402; -3.8185738483809524e9 -1.7123166484375e6 … -6.2063795222222224e7 119336.55902777778]), NFFTPlan with 284 sampling points for an input array of size(32, 32) and an output array of size(284,) with dims 1:2, ComplexF64[0.9304270810105981 - 0.1755525366839422im, 0.9999999857224657 - 6.37119066797087e-5im, 0.6968566352375541 - 0.5221791050707502im, 0.9304860627318433 + 0.17550202040945986im, 0.8028821861479044 + 0.3950138433443315im, 0.696909133455242 + 0.522149977275377im, 0.9999999859666728 - 6.454986502401626e-5im, 0.690100908057055 - 0.5259755617707915im, 0.9317353097306306 - 0.1646343527093735im, 0.9317937328057632 + 0.1645830409592267im … 0.8346090489457091 + 0.1795789988845693im, 0.834610800886739 + 0.17958571977595708im, 0.9181402138641586 + 0.09206131734480663im, 0.9120864418553256 + 0.2942607208902332im, 0.9120838428437319 + 0.29426742718619836im, 0.9559925804493471 + 0.09847981391058878im, 0.8881131676530709 + 0.289296040490046im, 0.9924358388500316 + 0.0028007013218525043im, 0.8881114735852544 + 0.28930284794286826im, 0.9999999999233905 + 7.3999334332420275e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]), gcache = JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9 … 274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP … :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256 … 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP) … (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25)), gcachep = JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22 … 275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im … 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9 … 274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im … 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP … :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379 … 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP) … (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)) … SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))]))We will also fix the total flux to be the observed value 1.1. This is because total flux is degenerate with a global shift in the gain amplitudes making the problem degenerate. To fix this we use the observed total flux as our value.
Moving onto our prior, we first focus on the instrument model priors. Each station requires its own prior on both the amplitudes and phases. For the amplitudes we assume that the gains are apriori well calibrated around unit gains (or 0 log gain amplitudes) which corresponds to no instrument corruption. The gain dispersion is then set to 10% for all stations except LMT, representing that we expect 10% deviations from scan-to-scan. For LMT we let the prior expand to 100% due to the known pointing issues LMT had in 2017.
using Distributions
using DistributionsAD
distamp = station_tuple(dvis, Normal(0.0, 0.1); LM = Normal(1.0))(AA = Distributions.Normal{Float64}(μ=0.0, σ=0.1), AP = Distributions.Normal{Float64}(μ=0.0, σ=0.1), AZ = Distributions.Normal{Float64}(μ=0.0, σ=0.1), JC = Distributions.Normal{Float64}(μ=0.0, σ=0.1), LM = Distributions.Normal{Float64}(μ=1.0, σ=1.0), PV = Distributions.Normal{Float64}(μ=0.0, σ=0.1), SM = Distributions.Normal{Float64}(μ=0.0, σ=0.1))For the phases, as mentioned above, we will use a segmented gain prior. This means that rather than the parameters being directly the gains, we fit the first gain for each site, and then the other parameters are the segmented gains compared to the previous time. To model this we break the gain phase prior into two parts. The first is the prior for the first observing timestamp of each site, distphase0, and the second is the prior for segmented gain ϵₜ from time i to i+1, given by distphase. For the EHT, we are dealing with pre-2*rand(rng, ndim) .- 1.5calibrated data, so often, the gain phase jumps from scan to scan are minor. As such, we can put a more informative prior on distphase.
We use AA (ALMA) as a reference station so we do not have to specify a gain prior for it.
distphase = station_tuple(dvis, DiagonalVonMises(0.0, inv(π^2)))(AA = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), AP = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), AZ = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), JC = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), LM = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), PV = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), SM = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688))In addition we want a reasonable guess for what the resolution of our image should be. For radio astronomy this is given by roughly the longest baseline in the image. To put this into pixel space we then divide by the pixel size.
beam = beamsize(dvis)
rat = (beam/(step(grid.X)))5.279832635689346To make the Gaussian Markov random field efficient we first precompute a bunch of quantities that allow us to scale things linearly with the number of image pixels. This drastically improves the usual N^3 scaling you get from usual Gaussian Processes.
crcache = MarkovRandomFieldCache(meanpr)VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}(sparse([1, 2, 32, 33, 993, 1, 2, 3, 34, 994 … 31, 991, 1022, 1023, 1024, 32, 992, 993, 1023, 1024], [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 58, 58, 58, 58, 58, 59, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64, 65, 65, 65, 65, 65, 66, 66, 66, 66, 66, 67, 67, 67, 67, 67, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 74, 74, 74, 74, 74, 75, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 77, 77, 77, 77, 78, 78, 78, 78, 78, 79, 79, 79, 79, 79, 80, 80, 80, 80, 80, 81, 81, 81, 81, 81, 82, 82, 82, 82, 82, 83, 83, 83, 83, 83, 84, 84, 84, 84, 84, 85, 85, 85, 85, 85, 86, 86, 86, 86, 86, 87, 87, 87, 87, 87, 88, 88, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 90, 90, 90, 91, 91, 91, 91, 91, 92, 92, 92, 92, 92, 93, 93, 93, 93, 93, 94, 94, 94, 94, 94, 95, 95, 95, 95, 95, 96, 96, 96, 96, 96, 97, 97, 97, 97, 97, 98, 98, 98, 98, 98, 99, 99, 99, 99, 99, 100, 100, 100, 100, 100, 101, 101, 101, 101, 101, 102, 102, 102, 102, 102, 103, 103, 103, 103, 103, 104, 104, 104, 104, 104, 105, 105, 105, 105, 105, 106, 106, 106, 106, 106, 107, 107, 107, 107, 107, 108, 108, 108, 108, 108, 109, 109, 109, 109, 109, 110, 110, 110, 110, 110, 111, 111, 111, 111, 111, 112, 112, 112, 112, 112, 113, 113, 113, 113, 113, 114, 114, 114, 114, 114, 115, 115, 115, 115, 115, 116, 116, 116, 116, 116, 117, 117, 117, 117, 117, 118, 118, 118, 118, 118, 119, 119, 119, 119, 119, 120, 120, 120, 120, 120, 121, 121, 121, 121, 121, 122, 122, 122, 122, 122, 123, 123, 123, 123, 123, 124, 124, 124, 124, 124, 125, 125, 125, 125, 125, 126, 126, 126, 126, 126, 127, 127, 127, 127, 127, 128, 128, 128, 128, 128, 129, 129, 129, 129, 129, 130, 130, 130, 130, 130, 131, 131, 131, 131, 131, 132, 132, 132, 132, 132, 133, 133, 133, 133, 133, 134, 134, 134, 134, 134, 135, 135, 135, 135, 135, 136, 136, 136, 136, 136, 137, 137, 137, 137, 137, 138, 138, 138, 138, 138, 139, 139, 139, 139, 139, 140, 140, 140, 140, 140, 141, 141, 141, 141, 141, 142, 142, 142, 142, 142, 143, 143, 143, 143, 143, 144, 144, 144, 144, 144, 145, 145, 145, 145, 145, 146, 146, 146, 146, 146, 147, 147, 147, 147, 147, 148, 148, 148, 148, 148, 149, 149, 149, 149, 149, 150, 150, 150, 150, 150, 151, 151, 151, 151, 151, 152, 152, 152, 152, 152, 153, 153, 153, 153, 153, 154, 154, 154, 154, 154, 155, 155, 155, 155, 155, 156, 156, 156, 156, 156, 157, 157, 157, 157, 157, 158, 158, 158, 158, 158, 159, 159, 159, 159, 159, 160, 160, 160, 160, 160, 161, 161, 161, 161, 161, 162, 162, 162, 162, 162, 163, 163, 163, 163, 163, 164, 164, 164, 164, 164, 165, 165, 165, 165, 165, 166, 166, 166, 166, 166, 167, 167, 167, 167, 167, 168, 168, 168, 168, 168, 169, 169, 169, 169, 169, 170, 170, 170, 170, 170, 171, 171, 171, 171, 171, 172, 172, 172, 172, 172, 173, 173, 173, 173, 173, 174, 174, 174, 174, 174, 175, 175, 175, 175, 175, 176, 176, 176, 176, 176, 177, 177, 177, 177, 177, 178, 178, 178, 178, 178, 179, 179, 179, 179, 179, 180, 180, 180, 180, 180, 181, 181, 181, 181, 181, 182, 182, 182, 182, 182, 183, 183, 183, 183, 183, 184, 184, 184, 184, 184, 185, 185, 185, 185, 185, 186, 186, 186, 186, 186, 187, 187, 187, 187, 187, 188, 188, 188, 188, 188, 189, 189, 189, 189, 189, 190, 190, 190, 190, 190, 191, 191, 191, 191, 191, 192, 192, 192, 192, 192, 193, 193, 193, 193, 193, 194, 194, 194, 194, 194, 195, 195, 195, 195, 195, 196, 196, 196, 196, 196, 197, 197, 197, 197, 197, 198, 198, 198, 198, 198, 199, 199, 199, 199, 199, 200, 200, 200, 200, 200, 201, 201, 201, 201, 201, 202, 202, 202, 202, 202, 203, 203, 203, 203, 203, 204, 204, 204, 204, 204, 205, 205, 205, 205, 205, 206, 206, 206, 206, 206, 207, 207, 207, 207, 207, 208, 208, 208, 208, 208, 209, 209, 209, 209, 209, 210, 210, 210, 210, 210, 211, 211, 211, 211, 211, 212, 212, 212, 212, 212, 213, 213, 213, 213, 213, 214, 214, 214, 214, 214, 215, 215, 215, 215, 215, 216, 216, 216, 216, 216, 217, 217, 217, 217, 217, 218, 218, 218, 218, 218, 219, 219, 219, 219, 219, 220, 220, 220, 220, 220, 221, 221, 221, 221, 221, 222, 222, 222, 222, 222, 223, 223, 223, 223, 223, 224, 224, 224, 224, 224, 225, 225, 225, 225, 225, 226, 226, 226, 226, 226, 227, 227, 227, 227, 227, 228, 228, 228, 228, 228, 229, 229, 229, 229, 229, 230, 230, 230, 230, 230, 231, 231, 231, 231, 231, 232, 232, 232, 232, 232, 233, 233, 233, 233, 233, 234, 234, 234, 234, 234, 235, 235, 235, 235, 235, 236, 236, 236, 236, 236, 237, 237, 237, 237, 237, 238, 238, 238, 238, 238, 239, 239, 239, 239, 239, 240, 240, 240, 240, 240, 241, 241, 241, 241, 241, 242, 242, 242, 242, 242, 243, 243, 243, 243, 243, 244, 244, 244, 244, 244, 245, 245, 245, 245, 245, 246, 246, 246, 246, 246, 247, 247, 247, 247, 247, 248, 248, 248, 248, 248, 249, 249, 249, 249, 249, 250, 250, 250, 250, 250, 251, 251, 251, 251, 251, 252, 252, 252, 252, 252, 253, 253, 253, 253, 253, 254, 254, 254, 254, 254, 255, 255, 255, 255, 255, 256, 256, 256, 256, 256, 257, 257, 257, 257, 257, 258, 258, 258, 258, 258, 259, 259, 259, 259, 259, 260, 260, 260, 260, 260, 261, 261, 261, 261, 261, 262, 262, 262, 262, 262, 263, 263, 263, 263, 263, 264, 264, 264, 264, 264, 265, 265, 265, 265, 265, 266, 266, 266, 266, 266, 267, 267, 267, 267, 267, 268, 268, 268, 268, 268, 269, 269, 269, 269, 269, 270, 270, 270, 270, 270, 271, 271, 271, 271, 271, 272, 272, 272, 272, 272, 273, 273, 273, 273, 273, 274, 274, 274, 274, 274, 275, 275, 275, 275, 275, 276, 276, 276, 276, 276, 277, 277, 277, 277, 277, 278, 278, 278, 278, 278, 279, 279, 279, 279, 279, 280, 280, 280, 280, 280, 281, 281, 281, 281, 281, 282, 282, 282, 282, 282, 283, 283, 283, 283, 283, 284, 284, 284, 284, 284, 285, 285, 285, 285, 285, 286, 286, 286, 286, 286, 287, 287, 287, 287, 287, 288, 288, 288, 288, 288, 289, 289, 289, 289, 289, 290, 290, 290, 290, 290, 291, 291, 291, 291, 291, 292, 292, 292, 292, 292, 293, 293, 293, 293, 293, 294, 294, 294, 294, 294, 295, 295, 295, 295, 295, 296, 296, 296, 296, 296, 297, 297, 297, 297, 297, 298, 298, 298, 298, 298, 299, 299, 299, 299, 299, 300, 300, 300, 300, 300, 301, 301, 301, 301, 301, 302, 302, 302, 302, 302, 303, 303, 303, 303, 303, 304, 304, 304, 304, 304, 305, 305, 305, 305, 305, 306, 306, 306, 306, 306, 307, 307, 307, 307, 307, 308, 308, 308, 308, 308, 309, 309, 309, 309, 309, 310, 310, 310, 310, 310, 311, 311, 311, 311, 311, 312, 312, 312, 312, 312, 313, 313, 313, 313, 313, 314, 314, 314, 314, 314, 315, 315, 315, 315, 315, 316, 316, 316, 316, 316, 317, 317, 317, 317, 317, 318, 318, 318, 318, 318, 319, 319, 319, 319, 319, 320, 320, 320, 320, 320, 321, 321, 321, 321, 321, 322, 322, 322, 322, 322, 323, 323, 323, 323, 323, 324, 324, 324, 324, 324, 325, 325, 325, 325, 325, 326, 326, 326, 326, 326, 327, 327, 327, 327, 327, 328, 328, 328, 328, 328, 329, 329, 329, 329, 329, 330, 330, 330, 330, 330, 331, 331, 331, 331, 331, 332, 332, 332, 332, 332, 333, 333, 333, 333, 333, 334, 334, 334, 334, 334, 335, 335, 335, 335, 335, 336, 336, 336, 336, 336, 337, 337, 337, 337, 337, 338, 338, 338, 338, 338, 339, 339, 339, 339, 339, 340, 340, 340, 340, 340, 341, 341, 341, 341, 341, 342, 342, 342, 342, 342, 343, 343, 343, 343, 343, 344, 344, 344, 344, 344, 345, 345, 345, 345, 345, 346, 346, 346, 346, 346, 347, 347, 347, 347, 347, 348, 348, 348, 348, 348, 349, 349, 349, 349, 349, 350, 350, 350, 350, 350, 351, 351, 351, 351, 351, 352, 352, 352, 352, 352, 353, 353, 353, 353, 353, 354, 354, 354, 354, 354, 355, 355, 355, 355, 355, 356, 356, 356, 356, 356, 357, 357, 357, 357, 357, 358, 358, 358, 358, 358, 359, 359, 359, 359, 359, 360, 360, 360, 360, 360, 361, 361, 361, 361, 361, 362, 362, 362, 362, 362, 363, 363, 363, 363, 363, 364, 364, 364, 364, 364, 365, 365, 365, 365, 365, 366, 366, 366, 366, 366, 367, 367, 367, 367, 367, 368, 368, 368, 368, 368, 369, 369, 369, 369, 369, 370, 370, 370, 370, 370, 371, 371, 371, 371, 371, 372, 372, 372, 372, 372, 373, 373, 373, 373, 373, 374, 374, 374, 374, 374, 375, 375, 375, 375, 375, 376, 376, 376, 376, 376, 377, 377, 377, 377, 377, 378, 378, 378, 378, 378, 379, 379, 379, 379, 379, 380, 380, 380, 380, 380, 381, 381, 381, 381, 381, 382, 382, 382, 382, 382, 383, 383, 383, 383, 383, 384, 384, 384, 384, 384, 385, 385, 385, 385, 385, 386, 386, 386, 386, 386, 387, 387, 387, 387, 387, 388, 388, 388, 388, 388, 389, 389, 389, 389, 389, 390, 390, 390, 390, 390, 391, 391, 391, 391, 391, 392, 392, 392, 392, 392, 393, 393, 393, 393, 393, 394, 394, 394, 394, 394, 395, 395, 395, 395, 395, 396, 396, 396, 396, 396, 397, 397, 397, 397, 397, 398, 398, 398, 398, 398, 399, 399, 399, 399, 399, 400, 400, 400, 400, 400, 401, 401, 401, 401, 401, 402, 402, 402, 402, 402, 403, 403, 403, 403, 403, 404, 404, 404, 404, 404, 405, 405, 405, 405, 405, 406, 406, 406, 406, 406, 407, 407, 407, 407, 407, 408, 408, 408, 408, 408, 409, 409, 409, 409, 409, 410, 410, 410, 410, 410, 411, 411, 411, 411, 411, 412, 412, 412, 412, 412, 413, 413, 413, 413, 413, 414, 414, 414, 414, 414, 415, 415, 415, 415, 415, 416, 416, 416, 416, 416, 417, 417, 417, 417, 417, 418, 418, 418, 418, 418, 419, 419, 419, 419, 419, 420, 420, 420, 420, 420, 421, 421, 421, 421, 421, 422, 422, 422, 422, 422, 423, 423, 423, 423, 423, 424, 424, 424, 424, 424, 425, 425, 425, 425, 425, 426, 426, 426, 426, 426, 427, 427, 427, 427, 427, 428, 428, 428, 428, 428, 429, 429, 429, 429, 429, 430, 430, 430, 430, 430, 431, 431, 431, 431, 431, 432, 432, 432, 432, 432, 433, 433, 433, 433, 433, 434, 434, 434, 434, 434, 435, 435, 435, 435, 435, 436, 436, 436, 436, 436, 437, 437, 437, 437, 437, 438, 438, 438, 438, 438, 439, 439, 439, 439, 439, 440, 440, 440, 440, 440, 441, 441, 441, 441, 441, 442, 442, 442, 442, 442, 443, 443, 443, 443, 443, 444, 444, 444, 444, 444, 445, 445, 445, 445, 445, 446, 446, 446, 446, 446, 447, 447, 447, 447, 447, 448, 448, 448, 448, 448, 449, 449, 449, 449, 449, 450, 450, 450, 450, 450, 451, 451, 451, 451, 451, 452, 452, 452, 452, 452, 453, 453, 453, 453, 453, 454, 454, 454, 454, 454, 455, 455, 455, 455, 455, 456, 456, 456, 456, 456, 457, 457, 457, 457, 457, 458, 458, 458, 458, 458, 459, 459, 459, 459, 459, 460, 460, 460, 460, 460, 461, 461, 461, 461, 461, 462, 462, 462, 462, 462, 463, 463, 463, 463, 463, 464, 464, 464, 464, 464, 465, 465, 465, 465, 465, 466, 466, 466, 466, 466, 467, 467, 467, 467, 467, 468, 468, 468, 468, 468, 469, 469, 469, 469, 469, 470, 470, 470, 470, 470, 471, 471, 471, 471, 471, 472, 472, 472, 472, 472, 473, 473, 473, 473, 473, 474, 474, 474, 474, 474, 475, 475, 475, 475, 475, 476, 476, 476, 476, 476, 477, 477, 477, 477, 477, 478, 478, 478, 478, 478, 479, 479, 479, 479, 479, 480, 480, 480, 480, 480, 481, 481, 481, 481, 481, 482, 482, 482, 482, 482, 483, 483, 483, 483, 483, 484, 484, 484, 484, 484, 485, 485, 485, 485, 485, 486, 486, 486, 486, 486, 487, 487, 487, 487, 487, 488, 488, 488, 488, 488, 489, 489, 489, 489, 489, 490, 490, 490, 490, 490, 491, 491, 491, 491, 491, 492, 492, 492, 492, 492, 493, 493, 493, 493, 493, 494, 494, 494, 494, 494, 495, 495, 495, 495, 495, 496, 496, 496, 496, 496, 497, 497, 497, 497, 497, 498, 498, 498, 498, 498, 499, 499, 499, 499, 499, 500, 500, 500, 500, 500, 501, 501, 501, 501, 501, 502, 502, 502, 502, 502, 503, 503, 503, 503, 503, 504, 504, 504, 504, 504, 505, 505, 505, 505, 505, 506, 506, 506, 506, 506, 507, 507, 507, 507, 507, 508, 508, 508, 508, 508, 509, 509, 509, 509, 509, 510, 510, 510, 510, 510, 511, 511, 511, 511, 511, 512, 512, 512, 512, 512, 513, 513, 513, 513, 513, 514, 514, 514, 514, 514, 515, 515, 515, 515, 515, 516, 516, 516, 516, 516, 517, 517, 517, 517, 517, 518, 518, 518, 518, 518, 519, 519, 519, 519, 519, 520, 520, 520, 520, 520, 521, 521, 521, 521, 521, 522, 522, 522, 522, 522, 523, 523, 523, 523, 523, 524, 524, 524, 524, 524, 525, 525, 525, 525, 525, 526, 526, 526, 526, 526, 527, 527, 527, 527, 527, 528, 528, 528, 528, 528, 529, 529, 529, 529, 529, 530, 530, 530, 530, 530, 531, 531, 531, 531, 531, 532, 532, 532, 532, 532, 533, 533, 533, 533, 533, 534, 534, 534, 534, 534, 535, 535, 535, 535, 535, 536, 536, 536, 536, 536, 537, 537, 537, 537, 537, 538, 538, 538, 538, 538, 539, 539, 539, 539, 539, 540, 540, 540, 540, 540, 541, 541, 541, 541, 541, 542, 542, 542, 542, 542, 543, 543, 543, 543, 543, 544, 544, 544, 544, 544, 545, 545, 545, 545, 545, 546, 546, 546, 546, 546, 547, 547, 547, 547, 547, 548, 548, 548, 548, 548, 549, 549, 549, 549, 549, 550, 550, 550, 550, 550, 551, 551, 551, 551, 551, 552, 552, 552, 552, 552, 553, 553, 553, 553, 553, 554, 554, 554, 554, 554, 555, 555, 555, 555, 555, 556, 556, 556, 556, 556, 557, 557, 557, 557, 557, 558, 558, 558, 558, 558, 559, 559, 559, 559, 559, 560, 560, 560, 560, 560, 561, 561, 561, 561, 561, 562, 562, 562, 562, 562, 563, 563, 563, 563, 563, 564, 564, 564, 564, 564, 565, 565, 565, 565, 565, 566, 566, 566, 566, 566, 567, 567, 567, 567, 567, 568, 568, 568, 568, 568, 569, 569, 569, 569, 569, 570, 570, 570, 570, 570, 571, 571, 571, 571, 571, 572, 572, 572, 572, 572, 573, 573, 573, 573, 573, 574, 574, 574, 574, 574, 575, 575, 575, 575, 575, 576, 576, 576, 576, 576, 577, 577, 577, 577, 577, 578, 578, 578, 578, 578, 579, 579, 579, 579, 579, 580, 580, 580, 580, 580, 581, 581, 581, 581, 581, 582, 582, 582, 582, 582, 583, 583, 583, 583, 583, 584, 584, 584, 584, 584, 585, 585, 585, 585, 585, 586, 586, 586, 586, 586, 587, 587, 587, 587, 587, 588, 588, 588, 588, 588, 589, 589, 589, 589, 589, 590, 590, 590, 590, 590, 591, 591, 591, 591, 591, 592, 592, 592, 592, 592, 593, 593, 593, 593, 593, 594, 594, 594, 594, 594, 595, 595, 595, 595, 595, 596, 596, 596, 596, 596, 597, 597, 597, 597, 597, 598, 598, 598, 598, 598, 599, 599, 599, 599, 599, 600, 600, 600, 600, 600, 601, 601, 601, 601, 601, 602, 602, 602, 602, 602, 603, 603, 603, 603, 603, 604, 604, 604, 604, 604, 605, 605, 605, 605, 605, 606, 606, 606, 606, 606, 607, 607, 607, 607, 607, 608, 608, 608, 608, 608, 609, 609, 609, 609, 609, 610, 610, 610, 610, 610, 611, 611, 611, 611, 611, 612, 612, 612, 612, 612, 613, 613, 613, 613, 613, 614, 614, 614, 614, 614, 615, 615, 615, 615, 615, 616, 616, 616, 616, 616, 617, 617, 617, 617, 617, 618, 618, 618, 618, 618, 619, 619, 619, 619, 619, 620, 620, 620, 620, 620, 621, 621, 621, 621, 621, 622, 622, 622, 622, 622, 623, 623, 623, 623, 623, 624, 624, 624, 624, 624, 625, 625, 625, 625, 625, 626, 626, 626, 626, 626, 627, 627, 627, 627, 627, 628, 628, 628, 628, 628, 629, 629, 629, 629, 629, 630, 630, 630, 630, 630, 631, 631, 631, 631, 631, 632, 632, 632, 632, 632, 633, 633, 633, 633, 633, 634, 634, 634, 634, 634, 635, 635, 635, 635, 635, 636, 636, 636, 636, 636, 637, 637, 637, 637, 637, 638, 638, 638, 638, 638, 639, 639, 639, 639, 639, 640, 640, 640, 640, 640, 641, 641, 641, 641, 641, 642, 642, 642, 642, 642, 643, 643, 643, 643, 643, 644, 644, 644, 644, 644, 645, 645, 645, 645, 645, 646, 646, 646, 646, 646, 647, 647, 647, 647, 647, 648, 648, 648, 648, 648, 649, 649, 649, 649, 649, 650, 650, 650, 650, 650, 651, 651, 651, 651, 651, 652, 652, 652, 652, 652, 653, 653, 653, 653, 653, 654, 654, 654, 654, 654, 655, 655, 655, 655, 655, 656, 656, 656, 656, 656, 657, 657, 657, 657, 657, 658, 658, 658, 658, 658, 659, 659, 659, 659, 659, 660, 660, 660, 660, 660, 661, 661, 661, 661, 661, 662, 662, 662, 662, 662, 663, 663, 663, 663, 663, 664, 664, 664, 664, 664, 665, 665, 665, 665, 665, 666, 666, 666, 666, 666, 667, 667, 667, 667, 667, 668, 668, 668, 668, 668, 669, 669, 669, 669, 669, 670, 670, 670, 670, 670, 671, 671, 671, 671, 671, 672, 672, 672, 672, 672, 673, 673, 673, 673, 673, 674, 674, 674, 674, 674, 675, 675, 675, 675, 675, 676, 676, 676, 676, 676, 677, 677, 677, 677, 677, 678, 678, 678, 678, 678, 679, 679, 679, 679, 679, 680, 680, 680, 680, 680, 681, 681, 681, 681, 681, 682, 682, 682, 682, 682, 683, 683, 683, 683, 683, 684, 684, 684, 684, 684, 685, 685, 685, 685, 685, 686, 686, 686, 686, 686, 687, 687, 687, 687, 687, 688, 688, 688, 688, 688, 689, 689, 689, 689, 689, 690, 690, 690, 690, 690, 691, 691, 691, 691, 691, 692, 692, 692, 692, 692, 693, 693, 693, 693, 693, 694, 694, 694, 694, 694, 695, 695, 695, 695, 695, 696, 696, 696, 696, 696, 697, 697, 697, 697, 697, 698, 698, 698, 698, 698, 699, 699, 699, 699, 699, 700, 700, 700, 700, 700, 701, 701, 701, 701, 701, 702, 702, 702, 702, 702, 703, 703, 703, 703, 703, 704, 704, 704, 704, 704, 705, 705, 705, 705, 705, 706, 706, 706, 706, 706, 707, 707, 707, 707, 707, 708, 708, 708, 708, 708, 709, 709, 709, 709, 709, 710, 710, 710, 710, 710, 711, 711, 711, 711, 711, 712, 712, 712, 712, 712, 713, 713, 713, 713, 713, 714, 714, 714, 714, 714, 715, 715, 715, 715, 715, 716, 716, 716, 716, 716, 717, 717, 717, 717, 717, 718, 718, 718, 718, 718, 719, 719, 719, 719, 719, 720, 720, 720, 720, 720, 721, 721, 721, 721, 721, 722, 722, 722, 722, 722, 723, 723, 723, 723, 723, 724, 724, 724, 724, 724, 725, 725, 725, 725, 725, 726, 726, 726, 726, 726, 727, 727, 727, 727, 727, 728, 728, 728, 728, 728, 729, 729, 729, 729, 729, 730, 730, 730, 730, 730, 731, 731, 731, 731, 731, 732, 732, 732, 732, 732, 733, 733, 733, 733, 733, 734, 734, 734, 734, 734, 735, 735, 735, 735, 735, 736, 736, 736, 736, 736, 737, 737, 737, 737, 737, 738, 738, 738, 738, 738, 739, 739, 739, 739, 739, 740, 740, 740, 740, 740, 741, 741, 741, 741, 741, 742, 742, 742, 742, 742, 743, 743, 743, 743, 743, 744, 744, 744, 744, 744, 745, 745, 745, 745, 745, 746, 746, 746, 746, 746, 747, 747, 747, 747, 747, 748, 748, 748, 748, 748, 749, 749, 749, 749, 749, 750, 750, 750, 750, 750, 751, 751, 751, 751, 751, 752, 752, 752, 752, 752, 753, 753, 753, 753, 753, 754, 754, 754, 754, 754, 755, 755, 755, 755, 755, 756, 756, 756, 756, 756, 757, 757, 757, 757, 757, 758, 758, 758, 758, 758, 759, 759, 759, 759, 759, 760, 760, 760, 760, 760, 761, 761, 761, 761, 761, 762, 762, 762, 762, 762, 763, 763, 763, 763, 763, 764, 764, 764, 764, 764, 765, 765, 765, 765, 765, 766, 766, 766, 766, 766, 767, 767, 767, 767, 767, 768, 768, 768, 768, 768, 769, 769, 769, 769, 769, 770, 770, 770, 770, 770, 771, 771, 771, 771, 771, 772, 772, 772, 772, 772, 773, 773, 773, 773, 773, 774, 774, 774, 774, 774, 775, 775, 775, 775, 775, 776, 776, 776, 776, 776, 777, 777, 777, 777, 777, 778, 778, 778, 778, 778, 779, 779, 779, 779, 779, 780, 780, 780, 780, 780, 781, 781, 781, 781, 781, 782, 782, 782, 782, 782, 783, 783, 783, 783, 783, 784, 784, 784, 784, 784, 785, 785, 785, 785, 785, 786, 786, 786, 786, 786, 787, 787, 787, 787, 787, 788, 788, 788, 788, 788, 789, 789, 789, 789, 789, 790, 790, 790, 790, 790, 791, 791, 791, 791, 791, 792, 792, 792, 792, 792, 793, 793, 793, 793, 793, 794, 794, 794, 794, 794, 795, 795, 795, 795, 795, 796, 796, 796, 796, 796, 797, 797, 797, 797, 797, 798, 798, 798, 798, 798, 799, 799, 799, 799, 799, 800, 800, 800, 800, 800, 801, 801, 801, 801, 801, 802, 802, 802, 802, 802, 803, 803, 803, 803, 803, 804, 804, 804, 804, 804, 805, 805, 805, 805, 805, 806, 806, 806, 806, 806, 807, 807, 807, 807, 807, 808, 808, 808, 808, 808, 809, 809, 809, 809, 809, 810, 810, 810, 810, 810, 811, 811, 811, 811, 811, 812, 812, 812, 812, 812, 813, 813, 813, 813, 813, 814, 814, 814, 814, 814, 815, 815, 815, 815, 815, 816, 816, 816, 816, 816, 817, 817, 817, 817, 817, 818, 818, 818, 818, 818, 819, 819, 819, 819, 819, 820, 820, 820, 820, 820, 821, 821, 821, 821, 821, 822, 822, 822, 822, 822, 823, 823, 823, 823, 823, 824, 824, 824, 824, 824, 825, 825, 825, 825, 825, 826, 826, 826, 826, 826, 827, 827, 827, 827, 827, 828, 828, 828, 828, 828, 829, 829, 829, 829, 829, 830, 830, 830, 830, 830, 831, 831, 831, 831, 831, 832, 832, 832, 832, 832, 833, 833, 833, 833, 833, 834, 834, 834, 834, 834, 835, 835, 835, 835, 835, 836, 836, 836, 836, 836, 837, 837, 837, 837, 837, 838, 838, 838, 838, 838, 839, 839, 839, 839, 839, 840, 840, 840, 840, 840, 841, 841, 841, 841, 841, 842, 842, 842, 842, 842, 843, 843, 843, 843, 843, 844, 844, 844, 844, 844, 845, 845, 845, 845, 845, 846, 846, 846, 846, 846, 847, 847, 847, 847, 847, 848, 848, 848, 848, 848, 849, 849, 849, 849, 849, 850, 850, 850, 850, 850, 851, 851, 851, 851, 851, 852, 852, 852, 852, 852, 853, 853, 853, 853, 853, 854, 854, 854, 854, 854, 855, 855, 855, 855, 855, 856, 856, 856, 856, 856, 857, 857, 857, 857, 857, 858, 858, 858, 858, 858, 859, 859, 859, 859, 859, 860, 860, 860, 860, 860, 861, 861, 861, 861, 861, 862, 862, 862, 862, 862, 863, 863, 863, 863, 863, 864, 864, 864, 864, 864, 865, 865, 865, 865, 865, 866, 866, 866, 866, 866, 867, 867, 867, 867, 867, 868, 868, 868, 868, 868, 869, 869, 869, 869, 869, 870, 870, 870, 870, 870, 871, 871, 871, 871, 871, 872, 872, 872, 872, 872, 873, 873, 873, 873, 873, 874, 874, 874, 874, 874, 875, 875, 875, 875, 875, 876, 876, 876, 876, 876, 877, 877, 877, 877, 877, 878, 878, 878, 878, 878, 879, 879, 879, 879, 879, 880, 880, 880, 880, 880, 881, 881, 881, 881, 881, 882, 882, 882, 882, 882, 883, 883, 883, 883, 883, 884, 884, 884, 884, 884, 885, 885, 885, 885, 885, 886, 886, 886, 886, 886, 887, 887, 887, 887, 887, 888, 888, 888, 888, 888, 889, 889, 889, 889, 889, 890, 890, 890, 890, 890, 891, 891, 891, 891, 891, 892, 892, 892, 892, 892, 893, 893, 893, 893, 893, 894, 894, 894, 894, 894, 895, 895, 895, 895, 895, 896, 896, 896, 896, 896, 897, 897, 897, 897, 897, 898, 898, 898, 898, 898, 899, 899, 899, 899, 899, 900, 900, 900, 900, 900, 901, 901, 901, 901, 901, 902, 902, 902, 902, 902, 903, 903, 903, 903, 903, 904, 904, 904, 904, 904, 905, 905, 905, 905, 905, 906, 906, 906, 906, 906, 907, 907, 907, 907, 907, 908, 908, 908, 908, 908, 909, 909, 909, 909, 909, 910, 910, 910, 910, 910, 911, 911, 911, 911, 911, 912, 912, 912, 912, 912, 913, 913, 913, 913, 913, 914, 914, 914, 914, 914, 915, 915, 915, 915, 915, 916, 916, 916, 916, 916, 917, 917, 917, 917, 917, 918, 918, 918, 918, 918, 919, 919, 919, 919, 919, 920, 920, 920, 920, 920, 921, 921, 921, 921, 921, 922, 922, 922, 922, 922, 923, 923, 923, 923, 923, 924, 924, 924, 924, 924, 925, 925, 925, 925, 925, 926, 926, 926, 926, 926, 927, 927, 927, 927, 927, 928, 928, 928, 928, 928, 929, 929, 929, 929, 929, 930, 930, 930, 930, 930, 931, 931, 931, 931, 931, 932, 932, 932, 932, 932, 933, 933, 933, 933, 933, 934, 934, 934, 934, 934, 935, 935, 935, 935, 935, 936, 936, 936, 936, 936, 937, 937, 937, 937, 937, 938, 938, 938, 938, 938, 939, 939, 939, 939, 939, 940, 940, 940, 940, 940, 941, 941, 941, 941, 941, 942, 942, 942, 942, 942, 943, 943, 943, 943, 943, 944, 944, 944, 944, 944, 945, 945, 945, 945, 945, 946, 946, 946, 946, 946, 947, 947, 947, 947, 947, 948, 948, 948, 948, 948, 949, 949, 949, 949, 949, 950, 950, 950, 950, 950, 951, 951, 951, 951, 951, 952, 952, 952, 952, 952, 953, 953, 953, 953, 953, 954, 954, 954, 954, 954, 955, 955, 955, 955, 955, 956, 956, 956, 956, 956, 957, 957, 957, 957, 957, 958, 958, 958, 958, 958, 959, 959, 959, 959, 959, 960, 960, 960, 960, 960, 961, 961, 961, 961, 961, 962, 962, 962, 962, 962, 963, 963, 963, 963, 963, 964, 964, 964, 964, 964, 965, 965, 965, 965, 965, 966, 966, 966, 966, 966, 967, 967, 967, 967, 967, 968, 968, 968, 968, 968, 969, 969, 969, 969, 969, 970, 970, 970, 970, 970, 971, 971, 971, 971, 971, 972, 972, 972, 972, 972, 973, 973, 973, 973, 973, 974, 974, 974, 974, 974, 975, 975, 975, 975, 975, 976, 976, 976, 976, 976, 977, 977, 977, 977, 977, 978, 978, 978, 978, 978, 979, 979, 979, 979, 979, 980, 980, 980, 980, 980, 981, 981, 981, 981, 981, 982, 982, 982, 982, 982, 983, 983, 983, 983, 983, 984, 984, 984, 984, 984, 985, 985, 985, 985, 985, 986, 986, 986, 986, 986, 987, 987, 987, 987, 987, 988, 988, 988, 988, 988, 989, 989, 989, 989, 989, 990, 990, 990, 990, 990, 991, 991, 991, 991, 991, 992, 992, 992, 992, 992, 993, 993, 993, 993, 993, 994, 994, 994, 994, 994, 995, 995, 995, 995, 995, 996, 996, 996, 996, 996, 997, 997, 997, 997, 997, 998, 998, 998, 998, 998, 999, 999, 999, 999, 999, 1000, 1000, 1000, 1000, 1000, 1001, 1001, 1001, 1001, 1001, 1002, 1002, 1002, 1002, 1002, 1003, 1003, 1003, 1003, 1003, 1004, 1004, 1004, 1004, 1004, 1005, 1005, 1005, 1005, 1005, 1006, 1006, 1006, 1006, 1006, 1007, 1007, 1007, 1007, 1007, 1008, 1008, 1008, 1008, 1008, 1009, 1009, 1009, 1009, 1009, 1010, 1010, 1010, 1010, 1010, 1011, 1011, 1011, 1011, 1011, 1012, 1012, 1012, 1012, 1012, 1013, 1013, 1013, 1013, 1013, 1014, 1014, 1014, 1014, 1014, 1015, 1015, 1015, 1015, 1015, 1016, 1016, 1016, 1016, 1016, 1017, 1017, 1017, 1017, 1017, 1018, 1018, 1018, 1018, 1018, 1019, 1019, 1019, 1019, 1019, 1020, 1020, 1020, 1020, 1020, 1021, 1021, 1021, 1021, 1021, 1022, 1022, 1022, 1022, 1022, 1023, 1023, 1023, 1023, 1023, 1024, 1024, 1024, 1024, 1024], [4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0 … -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0], 1024, 1024), [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0], [0.0 0.03842943919353914 … 0.15224093497742697 0.03842943919353936; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785; … ; 0.15224093497742697 0.1906703741709661 … 0.30448186995485393 0.19067037417096633; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785])One of the benefits of the Bayesian approach is that we can fit for the hyperparameters of our prior/regularizers unlike traditional RML appraoches. To construct this heirarchical prior we will first make a map that takes in our regularizer hyperparameters and returns the image prior given those hyperparameters.
fmap = let crcache=crcache
x->GaussMarkovRandomField(x, 1.0, crcache)
end#1 (generic function with 1 method)Now we can finally form our image prior. For this we use a heirarchical prior where the inverse correlation length is given by a Half-Normal distribution whose peak is at zero and standard deviation is 0.1/rat where recall rat is the beam size per pixel. For the variance of the random field we use another half normal prior with standard deviation 0.1. The reason we use the half-normal priors is to prefer "simple" structures. Gaussian Markov random fields are extremly flexible models, and to prevent overfitting it is common to use priors that penalize complexity. Therefore, we want to use priors that enforce similarity to our mean image. If the data wants more complexity then it will drive us away from the prior.
cprior = HierarchicalPrior(fmap, InverseGamma(1.0, -log(0.01*rat)))VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3399885194615155)
θ: 2.9412757865584154
)
)
We can now form our model parameter priors. Like our other imaging examples, we use a Dirichlet prior for our image pixels. For the log gain amplitudes, we use the CalPrior which automatically constructs the prior for the given jones cache gcache.
prior = NamedDist(
fg = Uniform(0.0, 1.0),
σimg = truncated(Normal(0.0, 1.0); lower=0.01),
c = cprior,
lgamp = CalPrior(distamp, gcache),
gphase = CalPrior(distphase, gcachep),
)VLBIImagePriors.NamedDist{(:fg, :σimg, :c, :lgamp, :gphase), Tuple{Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}(
dists: (Distributions.Uniform{Float64}(a=0.0, b=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=1.0); lower=0.01), VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3399885194615155)
θ: 2.9412757865584154
)
)
, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}(
dists: DiagNormal(
dim: 129
μ: [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
Σ: [0.010000000000000002 0.0 … 0.0 0.0; 0.0 0.010000000000000002 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 0.010000000000000002]
)
jcache: JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9 … 274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP … :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256 … 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP) … (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25))
)
, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}(
dists: VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}(
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
κ: [0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778 … 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778]
lnorm: -180.86855628209554
)
jcache: JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22 … 275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im … 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9 … 274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im … 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP … :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379 … 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP) … (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)) … SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))])
)
)
)
Putting it all together we form our likelihood and posterior objects for optimization and sampling.
lklhd = RadioLikelihood(sky, instrument, dvis; skymeta=metadata, instrumentmeta=metadata)
post = Posterior(lklhd, prior)Posterior{RadioLikelihood{Comrade.ModelMetadata{typeof(Main.sky), NamedTuple{(:ftot, :K, :meanpr, :grid, :cache, :gcache, :gcachep), Tuple{Float64, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, Matrix{Float64}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}, Comrade.ModelMetadata{typeof(Main.instrument), NamedTuple{(:ftot, :K, :meanpr, :grid, :cache, :gcache, :gcachep), Tuple{Float64, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, Matrix{Float64}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}, Tuple{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}}, Tuple{Comrade.ConditionedLikelihood{Comrade.var"#26#27"{Float64, Vector{Float64}}, Vector{ComplexF64}}}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, NamedTuple{(:U, :V, :T, :F), NTuple{4, Vector{Float64}}}}, VLBIImagePriors.NamedDist{(:fg, :σimg, :c, :lgamp, :gphase), Tuple{Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}}(RadioLikelihood
Number of data products: 1
, VLBIImagePriors.NamedDist{(:fg, :σimg, :c, :lgamp, :gphase), Tuple{Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}(
dists: (Distributions.Uniform{Float64}(a=0.0, b=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=1.0); lower=0.01), VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3399885194615155)
θ: 2.9412757865584154
)
)
, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}(
dists: DiagNormal(
dim: 129
μ: [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
Σ: [0.010000000000000002 0.0 … 0.0 0.0; 0.0 0.010000000000000002 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 0.010000000000000002]
)
jcache: JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9 … 274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP … :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256 … 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP) … (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25))
)
, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}(
dists: VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}(
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
κ: [0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778 … 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778]
lnorm: -180.86855628209554
)
jcache: JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22 … 275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im … 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9 … 274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im … 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP … :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379 … 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP) … (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)) … SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))])
)
)
)
)Reconstructing the Image and Instrument Effects
To sample from this posterior, it is convenient to move from our constrained parameter space to an unconstrained one (i.e., the support of the transformed posterior is (-∞, ∞)). This is done using the asflat function.
tpost = asflat(post)
ndim = dimension(tpost)1364Our Posterior and TransformedPosterior objects satisfy the LogDensityProblems interface. This allows us to easily switch between different AD backends and many of Julia's statistical inference packages use this interface as well.
using LogDensityProblemsAD
using Zygote
gtpost = ADgradient(Val(:Zygote), tpost)
x0 = randn(rng, ndim)
LogDensityProblemsAD.logdensity_and_gradient(gtpost, x0)(-3.703184201226103e7, [5.402397410750556e6, -2.7798738872739602e7, 8.783477104155637, 23.57396054555629, 7.883868445240969, 32.55261519903449, 8.240025410281866, 16.331219645584355, 849.8705665625981, 396.9041870545225 … 5925.529873295528, 3644.0042791216465, -7508.628060914201, -8265.628740423648, 22095.303838375075, 27246.98068687621, -199991.43415044196, -67093.98279354392, 23587.11200719565, -19385.479418322386])We can now also find the dimension of our posterior or the number of parameters we are going to sample.
This can often be different from what you would expect. This is especially true when using angular variables where we often artificially increase the dimension of the parameter space to make sampling easier.
To initialize our sampler we will use optimize using LBFGS
using ComradeOptimization
using OptimizationOptimJL
f = OptimizationFunction(tpost, Optimization.AutoZygote())
prob = Optimization.OptimizationProblem(f, prior_sample(rng, tpost), nothing)
ℓ = logdensityof(tpost)
sol = solve(prob, LBFGS(), maxiters=1_000, g_tol=1e-1);Now transform back to parameter space
xopt = transform(tpost, sol.u)(fg = 0.999999923585566, σimg = 0.010012394758643886, c = (params = [3.007056559992978e-7 -1.0767207514732477e-6 … -1.3360830896482375e-6 3.796015296104413e-7; 8.504438981335337e-7 5.919721634754177e-7 … 2.3476008985780363e-6 -7.708204448797021e-7; … ; -2.8091700250618577e-6 1.009444785036751e-6 … -3.0738163175338932e-6 4.014961187072578e-6; -5.164059585070686e-7 -2.3842794858603416e-7 … 1.686633262263826e-6 -1.7259378425010108e-6], hyperparams = 0.002897576471773542), lgamp = [-0.00010713838404966126, 0.0001511634371070609, 1.1288636925273015, -3.671408852413435e-5, 0.00030277487261706987, -0.0002504536910039823, 1.0967727997323906, -0.00016706603407498998, 0.00032095952005297405, -0.00014447950333392286 … 0.0004445392437167962, -0.0029709381633020904, 0.9762974174550304, -0.0029251788374876964, -0.00022406410267030263, 0.00029713883840577586, 2.814741083309631e-5, -0.0019237926216802742, 0.9089834782493547, -0.0013607216972725744], gphase = [-0.8136896808922691, 0.8975151341396079, -1.2440491949930408, -0.8704853129781674, -2.9034279710412108, -3.1392911815922058, -0.928181043712945, 1.3311329960599876, 0.9725646135734101, -0.9765841283485222 … -0.5844499358842421, 0.8340342085360564, -1.0301624088379864, 2.1131323309258176, 0.7643645503516703, -0.5241395521839227, -3.0764611223478706, 1.7119046831466715, -0.04262454210161644, -2.7625189740712766])Fitting gains tends to be very difficult, meaning that optimization can take a lot longer. The upside is that we usually get nicer images.
First we will evaluate our fit by plotting the residuals
using Plots
residual(vlbimodel(post, xopt), dvis)These look reasonable, although there may be some minor overfitting. This could be improved in a few ways, but that is beyond the goal of this quick tutorial. Plotting the image, we see that we have a much cleaner version of the closure-only image from Imaging a Black Hole using only Closure Quantities.
import WGLMakie as CM
img = intensitymap(skymodel(post, xopt), fovx, fovy, 128, 128)
CM.image(img, axis=(xreversed=true, aspect=1, title="MAP Image"), colormap=:afmhot)